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1  Introduction 


The  goal  of  this  research  was  to  study  the  computational  issues  involve  in  constructing  an  accurate,  efficient 
and  flexible  computational  framework  for  the  modeling  of  supersonics  reactive  flows,  and  to  apply  the  methods 
thus  developed  to  problems  originated  from  the  Air-Force  labs. 

The  main  components  of  this  framework  are: 

•  High  order  accuracy  numerical  methods  are  applied  in  order  to  resolve  the  important  detailed  features  of 
the  flow. 

-  Spectral  multidomain  methods. 

-  Weighted  Essentially  Non  Oscillatory  (WENO)  methods. 

-  High  Order  compact  finite  difference  schemes. 

•  Domain  decomposition  techniques  are  used, 

-  The  methods  are  based  on  penalty  impositions  of  interface  boundary  conditions  in  a  proven  time- 
stable  way.  The  discontinuous  Galerkin  methodology  is  implemented  in  the  finite  differences  compo¬ 
nent  of  the  research. 

•  Structured  as  well  as  unstructured  grid  multi-domain  spectral  methods  are  being  developed. 

-  In  the  structured  grid  methods  the  building  blocks  are  cubes  (in  3D)  or  rectangles  (in  2D).  The  grid 
is  predefined  by  a  tensor  product  of  Gaussian  points  to  assure  spectral  accuracy. 

-  The  standard  building  blocks  in  the  unstructured  grid  methods  are  the  tetrahedral  elements  (in  3D) 
or  triangles  (in  2D)  with  a  newly  developed  unstructured  nodal  basis  within  each  element. 

•  To  stabilize  the  spectral  and  the  compact  schemes: 

-  Superviscosity  terms  are  used  to  stabilize  the  methods  in  the  presence  of  shock  waves.  To  minimize 
the  amount  of  added  work,  the  superviscosity  method  is  carried  out  in  a  form  of  a  low  pass  exponential 
filter. 

-  Post-processing,  aimed  at  removing  the  Gibbs  phenomenon,  is  applied  to  recover  spectral  accuracy. 

•  Far  field  boundary  conditions  to  eliminate  non-physical  reflections  from  the  boundaries  of  the  computa¬ 
tional  domain,  are  applied. 

The  methods  thus  developed  are  being  currently  applied  to  the  simulations  of  recessed  cavity  flameholders 
and  to  the  study  of  mixing  and  instabilities  induced  by  shock  waves. 

In  the  following  we  will  review  past  accomplishments  and  future  tasks  to  be  performed.  It  composed  of  the 
following  Sections: 

•  Spectral  methods  for  discontinuous  problems. 

•  Recessed  cavity  flameholders 

•  Weighted  Essentially  Non-Oscillatory  (WENO)  schemes. 

•  Spectral  multi-domain  methods. 

•  Unstructured  grid  based  spectral  and  WENO  methods. 

•  Far  field  boundary  conditions  for  the  Euler  equations  of  Gas  Dynamics. 
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2  Spectral  Methods  for  Shock  Waves 

The  use  of  spectral  and  other  high  order  accuracy  methods  for  the  simulation  of  high  Mach  number  flows 
has  been  one  of  the  main  topics  of  this  research.  Spectral  methods,  due  to  their  accuracy,  resolve  small  scale 
phenomena  better  than  lower  order  methods  and  therefore  they  are  the  methods  of  choice  when  long  time 
integrations  are  performed  and  fine  details  of  the  flow  are  needed.  The  presence  of  large  gradients,  however, 
poses  difficulties  in  the  application  of  spectral  methods  both  on  the  theoretical  and  the  practical  level.  We 
would  like  to  describe  the  state  of  the  art  in  the  application  of  spectral  methods  to  discontinuous  problems. 
First,  we  will  briefly  review  the  main  components  of  the  spectral  shock  capturing  techniques,  and  will  list  the 
proposed  research  topics  in  this  area.  We  then  will  give  a  more  detailed  review  of  the  achievements  in  the  field. 
The  main  components  of  a  spectral  shock  capturing  techniques  are 

•  Spectral  viscosity  terms  are  added  to  stabilize  the  schemes.  The  following  methods  are  currently  being 
used: 


—  Spectral  viscosity,  where  the  viscosity  is  added  only  to  the  high  modes  of  the  equation. 

—  Spectral  super  viscosity.  Here  a  high  term  derivative  is  added  to  stabilize  the  equation  without 
changing  the  order. 

—  A  low  pass  filter,  equivalent  to  the  Spectral  Super  Viscosity  is  added. 

•  At  the  end  of  the  calculation  the  discontinuities  are  located.  There  are  several  methods  of  accomplishing 
that. 


-  Methods  based  on  decay  rate  of  the  coefficients. 

—  Methods  based  on  conjugate  sums. 

-  Wavelet  based  methods. 

•  Reconstruction  of  the  solution,  based  on  the  resolution  of  the  Gibbs  phenomenon  is  carried  out. 

There  are  many  practical  problems  in  the  above  outlined  framework  that  need  to  be  addressed. 

•  What  is  the  best  mechanism  for  stabilizing  the  spectral  scheme?  This  is  a  very  important  question  that 
has  a  far  reaching  impact  on  the  successful  simulations  of  reactive  flow.  We  propose  to  study  the  interplay 
between  stability  and  accuracy  and  to  find  the  optimal  way  of  stabilizing  the  spectral  scheme.  We  will 
study  the  spectral  stabilizers  in  the  context  of  real  life  applications. 

•  The  reconstruction  step  is  supported  by  a  solid  theoretical  foundation,  however  its  application  is  still  more 
an  art  than  science  in  the  sense  that  there  are  too  many  free  parameters  that  are  crucial  to  the  success  of 
the  method.  We  propose  to  study  the  role  of  those  free  parameters  in  order  to  come  with  automatic  post 
processing  step. 

The  application  of  spectral  methods  to  flows  containing  shock  waves  has  been  extensively  studied  by  us 
and  others  in  the  past  and  significant  advances  had  been  made.  We  briefly  summarize  the  current  state  of  the 
research. 

2.1  Approximation  Theory 

2.1.1  Enhanced  Convergence  by  the  Use  of  Filters 

One  manifestation  of  the  slow  and  nonuniform  convergence  of  the  spectral  approximation  X^u  of  a  piecewise 
smooth  function  u  is  found  in  the  linear  decay  of  the  global  expansion  coefficients,  This  realization  also 
suggests  that  one  could  attempt  to  modify  the  global  expansion  coefficients  to  enhance  the  convergence  rate  of 
the  spectral  approximation.  The  critical  question  to  address  naturally  becomes  exactly  how  one  should  modify 
the  expansion  to  ensure  enhanced  convergence  to  the  correct  solution. 

Let  us  consider  the  filtered  approximation.  Ti\ntjs;(x),  of  the  form 

N 

TnUk(x)=  ^  a  u„eyip{inx)  ,  (2.1) 

n=-N 
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where  Un  signifies  the  discrete  expansion  coefficients  of  Ujv(x.  t)  and  cr[r})  is  a  real  filter  function  with  the 
following  properties 


cr(-7?) 

£t(0)  =  1 

£r(9)(0)  =  0  l<9<2p-l 

.  cr(7?)  =0  |r?|  >  1 


(2.2) 


If  <7(7^)  has  at  least  2p  —  1  continuous  derivatives,  a(r})  is  termed  a  filter  of  order  2p. 

As  the  filter  is  nothing  more  than  a  low  pass  filter,  it  is  not  surprising  that  the  filtered  function  converges 
faster  than  the  unfiltered  original  expansion.  To  understand  exactly  how  much  the  filtering  modifies  the  con¬ 
vergence  rate,  let  us  assume  that  u{x)  is  piecewise  C'P[0,27r]  with  one  discontinuity  located  at  x  =  .$.  Let  us 
furthermore  assume  that  the  filter  is  of  order  2p.  Then  the  pointwise  error  of  the  filtered  approximation  is  given 
as 

|«(i)  -  , 

where  d{x,^)  measures  the  distance  from  x  to  the  point  of  discontinuity,  K{u)  is  uniformly  bounded  away 
from  the  discontinuity  and  a  function  of  u{x)  only.  Also  ||  •  ||L|[o,2ff]  signifies  the  broken  L2[0, 27r]-norm. 

While  the  details  of  the  proof  of  this  result  are  quite  technical  the  interpretation  of  the  result  is  simple, 
and  perhaps  somewhat  surprising.  It  states  the  convergence  rate  of  the  filtered  approximation  is  determined 
solely  by  the  the  order,  2p,  of  the  filter,  cr(p),  and  the  regularity  of  the  function,  u(x),  away  from  the  point 
of  discontinuity.  In  particular,  if  the  function,  w(x),  is  piecewise  analytic  and  the  order  of  the  filter  increases 
with  N,  one  recovers  an  exponentially  accurate  approximation  to  the  unfiltered  function  everywhere  except  very 
close  to  the  discontinuity. 

The  actual  choice  of  the  filter  function,  a{T]),  is  one  of  great  variety  and  numerous  alternatives  are  discussed 
in  the  literature.  A  particularly  popular  one  is  the  exponential  filter 


a{T])  =  exp  {-ar]^P)  , 

which  satisfies  all  the  conditions  in  Eq.(2.2)  except  that  of  being  zero  for  I77]  >  1.  However,  by  choosing 
a  =  -  Xhem,  with  em  representing  the  machine  accuracy,  a{Tj)  vanishes  for  all  practical  purposes  when  |7?1 
exceeds  one  and  the  exponential  filter  allows  for  the  recovery  of  a  piecewise  exponentially  accurate  approximation 
of  a  piecewise  analytic  function  away  from  the  point  of  discontinuity. 

When  applying  the  filter  it  is  worth  noticing  that,  as  for  differentiation,  the  spectral  filtering,  has  a  dual 
formulation,  expressed  as  a  matrix  operation,  of  the  form 

2N-\  1^1  n 

TNUNixi)  =  ^  FijU.^ixj)  ,  F,J  =  ^  X!  cos(n(xi  -  xj))  , 

j=0  n=0 

where  cq  =  cr^  ~  2  and  c„  =  I  otherwise. 

The  use  of  a  filter  in  a  Chebyshev  collocation  method  is  similar  to  that  of  the  Fourier  approximation,  i.e., 
it  can  be  expressed  as 

N 

J'ArUAr(x)  =  W„r„(x)  , 

n=0 

or  on  its  dual  form 


!F^U[ii{Xi)  —  ^  'F,jU;v(Xj)  ,  Fjj  =  —  'y  ^  (T  T„{Xi)Tn{Xj)  . 

j=0  ^ 

While  the  use  of  spectral  filtering  as  discussed  in  the  above  remains  the  most  popular  way  of  enhancing  the 
convergence  rate,  an  alternative  approach  can  be  realized  by  observing  that  the  Gibbs  phenomenon  is  also  a 
manifestation  of  the  global  nature  of  the  interpolating  polynomial.  In  other  words,  one  could  attempt  to  improve 
on  the  quality  of  the  approximation  by  localizing  the  approximation  close  to  the  point  of  the  discontinuity. 

This  approach,  known  as  physical  space  filtering,  operates  directly  on  the  interpolating  polynomials  rather 
than  the  expansion  coefficients.  To  illustrate  the  basic  idea  consider  the  filtered  Fourier  approximation 
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2A'-1 

Ti^Uixix)-  ^  u{Xj)lj{T)  , 
j=0 

where  xj  are  the  usual  equidistant  grid  and  the  filtered  Lagrange  interpolation  polynomial  takes  the  form 

=  ^  E  ^^(^)^exp(m(x-a:,))  . 

n=  —  N 

Clearly,  if  a{T])  —  1  we  recover  the  Fourier  interpolation  polynomial.  However,  we  do  not  need  to  maintain  the 
close  connection  to  the  trigonometric  interpolation  polynomials,  as  expressed  in  IJix),  but  can  choose  to  use 
any  reasonable  kernel,  tpiy.x),  that  approximates  a  delta  function  as  y  -  x  approaches  zero.  We  can  exemplify 
this  idea  by  considering 


1  ^ 

^{y,x)  =  pii{y,x))—  exp{in^{y,x))  , 


n^-N 


where 


■ 

We  assume  that  u{x)  e  —  £,?/  +  £]  and  p{i{y,x))  controls  the  amounts  of  localization  as 

P(0)  =  1  ,p(^(y,x))  =  0  for  1^(2/, x)|  >  1  , 

i.e.,  the  kernel  vanished  outside  of  the  symmetric  interval  [y  -  e,y  +  e].  Note  also  that  ^  YLn=-N 
is  an  approximation  to  a  y-centered  delta  function.  It  this  setting  it  had  been  shown  that  the  order  of  the  filter 
and  the  regularity  of  the  function  away  from  the  point  of  discontinuity  solely  determines  the  convergence  rate 
of  the  filtered  approximation.  Exponential  convergence  can  be  recovered  everywhere  except  very  close  to  the 
point  of  discontinuity  as  measured  through  e.  The  need  to  specify  the  size,  s,  of  the  symmetric  interval  remains 
a  practical  obstacle  to  the  use  of  the  physical  space  filters. 

While  the  use  of  filters  can  have  a  dramatic  impact  on  the  quality  of  the  convergence  of  a  global  approximation 
to  a  discontinuous  function,  such  techniques  are  unable  to  improve  on  the  quality  of  the  approximation  as  one 
approaches  the  point  of  discontinuity.  Moreover,  filtering  generally  treats  the  Gibbs  oscillations  as  a  source  of 
noise  and  attempts  to  remove  it.  However,  as  has  been  speculated  in  the  past  and  recently  shown  rigorously, 
the  Gibbs  oscillations  are  not  noise  but  rather  contain  sufficient  information  to  reconstruct  an  exponentially 
convergent  approximation  everywhere  provided  only  that  the  location  of  the  discontinuity  is  known,  i.e.,  the 
Gibbs  phenomenon  can  be  overcome  completely. 

2.1.2  Resolving  the  Gibbs  Phenomenon 

In  the  following  we  outline  the  key  elements  of  a  general  theory  that  establishes  the  possibility  of  recovering  a 
piecewise  exponentially  convergent  series  to  a  piecewise  analytic  function,  u{x)  G  L^,[— 1, 1],  having  knowledge 
of  the  global  expansion  coefficients  and  the  position  of  the  discontinuities  only. 

The  basic  element  of  this  new  approach  is  the  identification  of  a  new  basis  with  very  special  properties  and, 
subsequently,  the  expansion  of  the  slowly  convergent  truncated  global  expansion  in  this  new  basis.  Provided  this 
new  basis  satisfies  certain  conditions,  the  new  expansion  has  the  remarkable  property  that  it  is  exponentially 
convergent  to  the  original  piecewise  analytic  function  even  though  it  uses  information  from  the  slowly  convergent 
global  expansion. 

As  previously  we  assume  that 


N 

Vnu{x)  -  Y^Un(t>n{x) 

n=0 


Hri  —  (w,  (i>Ti)w 

In 


with  7n  =  Note  in  particular  that  this  also  contains  the  Fourier  case  provided 


o 


&n  {x)  -  exp 


Let  us  also  assume  that  there  exists  an  interval  [a,  6]  C  [—1,1]  in  which  u{x)  is  analytic  and,  furthermore, 
that  the  original  truncated  expansion  is  pointwise  convergent  in  all  of  [-1,1]  with  the  exception  of  a  finite 
number  of  points.  We  shall  introduce  the  scaled  variable 


ax)  =  -i+2^  . 

0  —  a 

Clearly,  ^  ;  [a,  6]  ->■  [-1, 1], 

We  define  a  new  basis,  which  is  orthogonal  in  the  weighted  inner  product,  (■,  )„.  where  A  signifies 

that  the  weight,  w{x),  may  depend  on  A,  i,e., 

—  Ilv^'nlli2,[-l,l]'^*n  —  ')ln^kn  ■ 

Furthermore,  we  require  that  if  v{^)  is  analytic  then 


'PMO  =  A 

n=0 

is  pointwise  exponentially  convergent  as  A  increases,  i.e., 


||u  — 'PAn||r,~[_],i]  <  Ce  , 

with  c  >  0. 

A  final  condition,  however,  sets  this  basis  apart  and  is  central  in  order  to  overcome  the  Gibbs  phenomenon. 
We  shall  require  that  there  exists  a  number  0  <  1,  such  that  for  A  =  0N  we  have 


(2,3) 


for  k  >  N,  n  <  X  and  q  <  1,  The  interpretation  of  this  condition  is  that  the  projection  of  the  high  modes  of 
(ph  onto  the  basis,  is  exponentially  small  in  the  interval,  ^  €  [—1, 1],  In  other  words,  by  reexpanding  the 
slowly  decaying  <;!l)i-based  global  expansion  in  the  local  ?^'^-basis  an  exponentially  accurate  local  approximation 
is  recovered.  Moreover,  this  can  be  achieved  everywhere  in  the  domain  where  u(x)  is  analytic.  We  shall  term 
this  latter,  and  crucial,  condition  on  tl’n  the  Gibbs  condition  to  emphasize  its  close  connection  to  the  resolution 
of  the  Gibbs  phenomenon. 

Provided  only  that  the  V'^-basis,  which  we  term  the  Gibbs  complementary  basis,  is  complete  we  recover  the 
key  result 


n=0 


<  C  exp(— cAf)  , 


where  A  =  0N  and  u(x)  is  analytic  in  the  interval  [a, 6]. 

In  other  words,  if  a  Gibbs  complementary  basis  exists  it  is  possible  to  reconstruct  a  piecewise  exponentially 
convergent  approximation  to  a  piecewise  analytic  function  from  the  information  contained  in  the  original  very 
slowly  converging  global  approximation  using  only  knowledge  about  the  location  of  the  points  of  discontinuity. 
Hence,  the  impact  of  the  Gibbs  phenomenon  can  be  overcome. 

A  constructive  approach  to  the  identification  of  the  complementary  basis  is  currently  unknown  and  subject 
to  future  research.  The  existence  of  such  a  basis,  however,  has  been  established  by  carefully  examining  the 
properties  of  the  basis 


where  represent  the  Gegenbauer  polynomials,  also  known  as  the  symmetric  .Jacobi  polynomials  or  the 

ultraspherical  polynomials.  We  proved  that  this  may  serve  as  a  complementary  Gibbs  basis  to  the  Fourier  baisis, 
the  Chebyshev  polynomials  and  the  Legendre  polynomials. 
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We  have  focused  on  the  resolution  of  the  Gibbs  phenomenon  for  continuous  spectral  expansions  which, 
unfortunately,  have  little  practical  importance.  A  similar  discussion  and  analysis  for  the  discrete  expansion  is 
complicated  by  the  introduction  of  the  aliasing  error.  The  proof  that  the  Gegenbauer  polynomials  remain  a 
Gibbs  complementary  basis  for  the  Fourier  and  Chebyshev  polynomials  has  be  also  completed 

2.2  Convergence  Results  for  Nonlinear  Hyperbolic  Problems 

The  return  to  the  general  nonlinear  problem 


du 


dfju) 

dx 


=  0 


(2.4) 


introduces  additional  issues  and  new  problems  which  are  not  relevant  or  of  less  importance  for  the  variable 
coefficient  problem. 

One  of  the  central  difficulties  in  using  spectral  methods  for  the  solution  of  nonlinear  conservation  laws 
lies  in  the  potential  development  of  nonsmooth  solutions  in  finite  time  even  for  problems  with  very  smooth 
initial  conditions.  As  we  have  discussed  previously,  this  introduces  the  Gibbs  phenomenon  which,  through  the 
nonlinearity,  interacts  with  the  solution.  What  we  shall  discuss  in  the  following  is  the  impact  this  has  on  the 
performance  of  the  numerical  approximation  and  techniques  that  allow  us  to  recover  accurate  and  physically 
meaningful  solutions  to  the  conservation  laws  even  when  the  Gibbs  oscillations  are  apparent. 


2.2.1  Stability  by  the  Use  of  Filters 


Maintaining  stability  of  the  numerical  approximation  becomes  increasingly  hard  as  the  discontinuity  evolves 
and  generates  energy  with  higher  and  higher  frequency  content.  This  process,  amplified  by  the  nonlinear  mixing 
of  the  Gibbs  oscillations  and  the  numerical  solution,  eventually  renders  the  scheme  unstable. 

Understanding  the  source  of  the  stability  problem,  i.e.,  accumulation  of  high  frequency  energy,  also  suggests 
a  possible  solution  by  introducing  a  dissipative  mechanism  that  continuously  remove  these  high  frequency 
components. 

A  classical  way  to  accomplish  this  is  to  modify  the  original  problem  by  adding  artificial  dissipation  as 


du  df{u) 
dt  ^  dx 


d^pu 

dx'^P 


A  direct  implementation  of  this,  however,  may  be  costly  and  could  introduce  additional  stiffness  which  would 
limit  the  stable  time  step.  We  shall  hence  seek  a  different  approach  to  reach  a  similar  goal. 

To  understand  the  impact  of  using  the  filter  at  regular  intervals  as  a  stabilizing  mechanism,  consider  the  use 
of  an  exponential  filter 


a{T))  —  exp  [-ar]-P)  . 

As  discussed  in  Sec.  2.1.1  this  filter  allows  for  a  dramatic  improvement  in  the  accuracy  of  the  approximation 
away  from  points  of  discontinuity. 

To  appreciate  its  impact  on  stability,  consider  the  generic  initial  value  problem 


with  the  pseudospectral  Fourier  approximation 


—  «  =  Cnu  . 
at 

Advancing  the  solution  from  t  =  0  to  t  =  At  followed  by  the  filtering  is  conveniently  expressed  as 

u(Af)  =  J'jv  exp(£;vAt)u(0)  . 

If  we  first  assume  that  represents  the  constant  coefficient  hyperbolic  problem,  i.e.,  £  =  a^,  we  recover  that 

u,.(A/)  =  exp(-a7j^'’ +  a{iA:)Af)u„(0)  ,  (2.5) 

i.e.,  we  are  in  fact  computing  the  solution  to  the  modified  problem 


du  dv  (~1)^  d'^u 
dt  ^  dx  ^  AtA^'^P  dx~P 

The  effect  of  the  filter  is  thus  equivalent  to  the  classical  approach  of  adding  a  small  dissipative  term  to 
the  original  equation,  but  the  process  of  adding  the  dissipation  is  very  simple.  Note  in  particular  that  AtN 
essentially  represents  the  CFL  condition  and  hence  is  of  order  one. 

For  a  general  £,  e.g.,  with  a  variable  coefficient  or  of  a  nonlinear  form,  for  which  !F,\  and  C,\  no  longer 
commute,  the  modified  equation  being  solved  takes  the  form 


du 

dt 


£u  —  a 


(-l)P  d^Pu 

AtN^d^ 


+  0{AP) 


by  viewing  the  application  of  the  filter  as  an  operator  splitting  problem  [16] 

With  this  in  mind  it  is  not  surprising  that  using  a  filter  has  a  stabilizing  effect.  Moreover,  we  observe 
that  if  p  increases  with  N  the  modification  caused  by  the  filter  vanishes  spectrally  as  N  increases.  These 
loose  arguments  for  the  stabilizing  effect  of  filtering  have  been  put  on  firm  ground  for  problem  with  smooth 
and  nonsmooth  initial  data  [77,  89]  for  the  Fourier  approximation  to  the  general  variable  coefficient  problem. 
These  results,  however,  are  typically  derived  under  the  assumption  that  a{ri)  is  of  polynomial  form.  While  such 
filtering  indeed  stabilizes  the  approximation  it  may  also  reduce  the  global  accuracy  of  the  scheme  [77,  42].  Let 
us  therefore  briefly  consider  the  stabilizing  effect  of  using  the  exponential  filter  in  the  pseudospectral  Fourier 
approximation  of  hyperbolic  equations,  known  to  be  only  weakly  unstable. 

Consider  the  filtered  approximation  on  the  form 


,  (2.6) 

where  the  superviscosity  term  on  the  right  can  be  implemented  through  a  filter  and 

To  establish  stability,  let  us  rewrite  Eq.(2.6)  as 


where 


dui\i 

dt 


Af\U!^  pf  Mzui^i  =  e/v(— 


d'^PUN 

dx'^P 


Af\Ut<  =  -—l!wa{x)uiM  +  -Iat  (0(1) 
is  the  skew-symmetric  form  of  the  operator, 

1 


du 


N 


dx 


and 


To  establish  stability,  consider 


i,r  ^-T  (  (  ^9un\  da{x)ujw 

=  -I, 


1  da{x)uN  1  3  „  ,  , 


2dt 


II^^wIIa'  =  -[wAf,MwAr]A'  -  [wn,-A/2Wa(]/V 

—  [uAf,  AsUAr]/)^  + 


UN,eN{-lf^^ 


d'^PuN 

dx'^P 


N 


Clearly  [ujw ,Ai\un]n  =  0  due  to  the  skew-symmetry  of  A\u,\’  and  by  inspection  we  can  bound 

[ua’,A/2W/v]a'  <  k  max  |aa:(x)|||uA’irA-  . 

■  3-€[0.27r] 
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It  is  indeed  the  term  associated  with  that  is  the  source  of  the  problems.  To  appreciate  this,  note  that 

if  Vn  was  used  rather  than  I/v  such  that  differentiation  and  truncation  commute,  the  term  would  vanishes 
identically  and  the  scheme  would  be  stable.  To  bound  this  term,  we  can  use  that 


[un,N'3Uim]n  <  C  (||u;v1|a  +  IpVsUivIl^)  . 

Noting  that  ||A/3UA|li2[o  27r]  nothing  more  than  the  commutation  error  and  that  ||  •  is  L^-equivalent,  we 
can  obtain 


[u;v,  A/snA'lAT  <  c  (||uN|li2|o,2x]  +  Ili2[0.27r] ) 

where  C  depends  on  a{x)  and  its  first  p  derivatives.  If  we  finally  note  that 


dx^p 


-eAllni^’llA  . 


J  N 


it  is  clear  that  we  can  always  choose  En  —  AN‘^~^p  and  A  sufficiently  large  to  ensure  stability.  In  other  words, 
using  an  exponential  filter  is  sufficient  to  stabilize  the  Fourier  approximation. 

There  is  one  central  difference  in  the  effect  of  using  the  filter  in  the  Fourier  and  the  Chebyshev  approximation. 
In  the  latter,  the  modified  equation  takes  the  form 


du 

dt 


£u  ~  a 


i-ir 

AtN^P 


u  +  OiAt^)  . 


(2.7) 


Hence,  while  the  filtering  continues  to  introduce  dissipation,  it  is  spatially  varying.  In  particular,  it  vanishes 
as  one  approaches  the  boundaries  of  the  domain.  In  computations  with  moving  discontinuities  this  may  be 
a  source  of  problems  since  the  stabilization  decreases  as  the  discontinuity  approaches  the  boundaries  of  the 
computational  domain. 


2.2.2  Spectrally  Vanishing  Viscosity  and  Entropy  Solutions 

The  foundation  of  a  convergence  theory  for  spectral  approximations  to  conservation  laws  has  been  laid  in 
[91,  75]  for  the  periodic  case  and  subsequently  extended  in  [76]  to  the  Legendre  approximation  and  recently  to 
the  Chebyshev-Legendre  scheme  in  [73,  74], 

To  appreciate  the  basic  elements  of  this  convergence  theory  let  us  first  restrict  ourselves  to  the  periodic  case. 
For  the  discrete  approximation  to  Eq.(2.4)  we  must  add  a  dissipative  term  that  is  strong  enough  to  stabilize 
the  approximation,  yet  small  enough  not  to  ruin  the  spectral  accuracy  of  the  scheme.  In  [91,  75]  the  following 
spectral  viscosity  method  was  considered 

where 

=  ^  exp{inx)  . 

m<|n|< 

To  ensure  that  stability  is  maintained  m  should  not  be  taken  too  big.  On  the  other  hand,  taking  m  too  small 
will  impact  the  accuracy  in 

a  negative  way.  An  acceptable  compromise  seems  to  be 

m  ~  ,  6>  <  . 

2p 

Moreover,  the  smoothing  factors,  should  only  be  activated  for  high  modes  as 


for  |n|  >  m  and  Q„  =  1  otherwise.  Finally,  we  shall  assume  that  the  amplitude  of  the  viscosity  is  small  as 


dxP 


Q m  £) 


dPuN 


dxP 


dPuN 

dx.P 


(2.8) 
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c 

~  jV2p-l  • 

Under  these  assumptions,  one  can  prove  for  p  =  1  that  the  solution  is  bounded  in  L°°[0,27r]  and  obtain  the 
estimate 


II^Af||L2[0,27r]  + 


duN 

dx 


<C  . 


Convergence  to  the  correct  entropy  solution  then  follows  from  compensated  compactness  arguments  [91,  75]. 

To  realize  the  close  connection  between  the  spectral  viscosity  method  and  the  use  of  filters  discussed  in  Sec. 
2.2.1,  consider  the  simple  case  where  f{u)  =  au.  In  this  case,  the  solution  to  Eq.(2.8)  is  given  as 


Unit)  =  exp  [inat  —  ENTi^Qnj  Un{0)  ,  |n|  >  m  , 

'which  is  equivalent  to  the  effect  of  the  filtering  discussed  in  Sec.  2.2.1.  Note  that  the  direct  application  of  the 
vanishing  viscosity  term  in  Eq.(2.8)  amounts  to  2p  spatial  derivatives  while  filtering  as  discussed  in  Sec.  2.1.1 
can  be  done  at  little  or  no  additional  cost. 

For  p  7^  1  a  bound  on  the  L°°[0,27r]  is  no  longer  known.  However,  experience  suggests  that  it  is  better  to 
filter  from  the  first  mode  but  to  employ  a  slower  decay  of  the  expansion  coefficients,  corresponding  to  taking 
p  >  1.  This  yields  the  superviscosity  method  in  which  one  solves 


dt 


d-T>u,^ 

dx'^p 


which  we  recognize  from  Eq.(2.6)  as  being  equivalent  to  that  obtained  when  using  a  high-order  exponential 
filter. 

The  vanishing  viscosity  approximation  to  Eq.(2.4)  using  a  Chebyshev  collocation  approach  takes  the  form 


where  again 


\/l  - 

ax 


2p 


UN  +  Bun 


£n  ~ 


C 


and  p  grows  with  N  [76].  Here  the  boundary  operator.  Bun,  may  vanish  or  it  may  take  values  as 


-(1  -x)T'^{x) 

27’;(-l) 


{uNi-l,t)  -  g  ) 


_+(l  +a:)r;,f(x) 


2rA,(i) 


(u;v(l,<)-p+) 


which  we  recognize  as  the  weakly  imposed  penalty  terms.  Note  again  that  the  vanishing  viscosity  term  is 
equivalent  to  that  obtained  from  the  analysis  of  the  effect  of  spectral  space  filtering,  Eq.(2.7).  Similar  results 
can  be  obtained  for  the  Legendre  approximation  and  for  the  Chebyshev-Legendre  method  for  which  convergence 
has  been  proven  [73,  74],  using  arguments  similar  to  those  in  [91,  75], 

using  arguments  similar  to  those  in  [91,  75],  for  p  =  1  as  well  as  for  p  >  1. 


2.2.3  Conservation 

It  is  natural  to  question  whether  the  introduction  of  the  artificial  Gibbs  oscillations  has  any  impact  on  the 
basic  physical  properties  described  by  the  conservation  law,  e.g.,  mass  conservation  and  the  speed  by  which 
discontinuities  propagate. 

To  come  to  an  understanding  of  this,  assume  a  spatially  periodic  problem  and  consider  the  pseudospectral 
Fourier  scheme 


dt 


u  -f  D/  =  0 
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where  u  =  [uA'(0,t).  ..,ua-(i2A'-i,<)]^  represent  the  grid  vector  and  the  interpolation  of  the  flux  is  given  as 
The  first  thing  to  note  is  that 

r2-  i'27r 

/  Ul\!{x,t)dx=  I  UN{x,0)dx  , 

Jo  Jo 

as  an  immediate  consequence  of  the  accuracy  of  the  trapezoidal  rule  and  the  assumption  of  periodicity.  Hence, 
the 

approximation  conserves  the  'mass’  of  the  interpolation  of  the  initial  conditions. 

Let  us  introduce  a  smooth  periodic  test  function,  tpix,  t),  with  the  corresponding  grid  vector,  tp  =  1),  ..,tp!\r{x2N~ 

The  test  function,  ip{x,t),  is  assumed  to  vanish  at  large  t.  If  we  consider  [35] 

^-(|„  +  D/)=0, 

and  utilize  the  accuracy  of  the  trapezoidal  rule  we  recover 


4)N{x,t) 


duN{x,t) 

dt 


dipN{x,t) 

dx 


1^Nf{uN{x,  t), 


t) 


dx  —  0  , 


after  integration  by  parts  which  is  permitted  if  the  solution,  ui\j{x,t),  is  bounded.  This  implicitly  assumes  that 
the  numerical  approximation  itself  is  stable  which  generally  implies  that  a  vanishing  viscosity  term  is  to  be 
added,  potentially  through  the  use  of  a  filter. 

Integrating  over  time  and  by  parts  once  more,  we  recover  the  result 


ui\r{x,t) 


dip!^{x,t) 

dt 


+ 


dxl}N{x,t) 

dx 


lNf{UN{x,t),t) 


dxdt 


f2n 

+  /  'tp;^{x,0)Ufj{x,0)  dx  =  0  . 

Jo 


Thus,  for  -t  00  the  solution,  uyv(x,t),  is  a  weak  solution  to  the  conservation  law.  This  essentially  implies  that 
the  limit  solution  satisfies  the  Rankine-Hugoniot  conditions  which  again  guarantees  that  shocks  propagate  at 
the  right  speed  to  within  the  order  of  the  scheme.  Results  similar  to  these  have  been  obtained  for  the  Chebyshev 
approximation  to  the  conservation  law  [35]. 

To  appreciate  that  the  addition  of  the  vanishing  viscosity  has  no  impact  on  the  conservation  of  the  scheme, 
consider  the  Legendre  superviscosity  case  [76]  and  let  be  a  test  function  in  C^[-l,  1]  that  vanishes  at 

the  endpoints.  Taking  V’n_i(x,<)  =  Iiw-\i!!{x,t),  then  clearly  ‘ip!\i-i{x,t)  — >  ■ip{x,t),  {tpN-i)x{x,t)  — t  tpx{x,t), 
and  {'tpN-\)t[x,t)  ‘ilii{x,t)  uniformly  in  N 

Since  ■0N-i(a:)  is  a  polynomial  that  vanishes  on  the  boundaries  we  have 


J  {l  +  x)P'^{x)'ipN^i{x,t)dx  =  0,  j  (l-x)F^, 


{x)ip{^'-.i{x,t)  dx  =  0. 


Moreover,  integration  by  parts  yields  that 

£n(-1)p 


lim 


.1)P 

—  J  xpjsj^i{x,t) 


dx  dx 


2p 


ui\!{x,  t)dx  =  0 


A^— +O0 

Hence,  the  superviscosity  term  does  not  cause  any  problems  and  one  can  show  that 

dxpN-i{x,t) 


^  ^  /  ,^d^pN-dx,t)  , 

UN(x,t) - -  + 


-a: 


dt 


lNf{UN{x,  t))- 


dx 


-j: 


dx  dt 


u/v(x,  0)j/j/v-i  (a;,  0)  dx  =  0  . 


The  main  conclusion  of  this  is  that  if  UA’(x,f)  is  a  solution  to  the  Legendre  collocation  approximation  at  the 
Gauss-Lobatto  points  and  if  ui\!{x.t)  converges  almost  everywhere  to  a  function  u{x,t),  then  u{x,t)  is  a  weak 
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solution  to  Eq.(2.4).  The  technical  details  of  this  proof  can  be  found  in  [14]  where  also  a  similar  result  for  the 
Chebyshev  superviscosity  approximation  is  given. 

The  theory  of  convergence  of  spectral  methods  equipped  with  spectral  viscosity  or  superviscosity  is  limited 
to  the  scalar  case  as  discussed  in  Sec.  2.2.2.  For  the  system  case  a  more  limited  result  can  be  obtained,  stating 
that  if  the  solution  converges  to  a  bounded  solution,  it  converges  to  the  correct  weak  solution. 
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3  Recessed  Cavity  Flameholders 

In  [29,  30]  we  applied  the  spectral  penalty  method  to  simulate  supersonic  combustion  problems  in  recessed 
cavities  in  order  to  establish  the  efficacy  of  recessed  cavity  flame-holders. 

We  considered  two  different  cases;  {!) Non- reactive  flows  with  two  chemical  species  and  {2)Reactive  flows 
with  four  chemical  species.  (3) Injector-cavity  flows  with  four  chemical  species. 

Recessed  cavities  provide  a  high  temperature,  low  speed  recirculating  region  that  can  support  the  production 
of  radicals  created  during  chemical  reactions.  This  stable  and  efficient  flame-holding  performance  by  the  cavity 
is  achieved  by  generating  a  recirculation  region  inside  the  cavity  where  a  hot  pool  of  radicals  forms  resulting  in 
reducing  the  induction  time  and  thus  obtaining  the  auto-ignition.  Experiments  have  shown  that  such  efficiency 
depends  on  the  geometry  of  the  cavity  such  as  the  degree  of  the  slantness  of  the  aft  wall  and  the  length  to 
depth  ratio  of  cavity  LfD.  Thus  one  can  optimize  the  flame-holding  performance  by  properly  adjusting  the 
geometrical  parameters  of  the  cavity  flame-holder  system  for  a  given  supersonic  flight  regime.  There  are  two 
major  issues  of  such  cavity  flame-holder  system  that  need  to  be  investigated:  {\.)What  is  the  optimal  angle  of 
the  aft  wall  for  a  given  L/D?  and  {2)How  does  the  fuel  injection  interact  with  cavity  flows?  An  answer  to  these 
questions  require  both  a  comprehensive  laboratory  and  numerical  experiments. 

In  this  work,  we  solved  the  full  compressible  Navier- Stokes  equations  with  chemical  reactions  without  any 
turbulence  model,  using  a  multi-domain  spectral  method. 

Results  of  several  numerical  studies  including  in  the  study  have  shown  that  the  stability  of  the  recirculation 
inside  cavity  is  enhanced  for  the  lower  angle  of  cavity  compared  to  the  rectangular  cavity.  The  present  study, 
however,  gives  more  accurate  and  finer  details  of  the  fields  than  those  done  by  lower  order  numerical  experiments. 
We  show  that  a  stationary  recirculation  region  is  not  formed  inside  the  cavity  contrary  to  what  the  lower  order 
schemes  predict.  A  quantitative  analysis  made  in  this  study  shows  that  the  lower  angled  wall  of  the  cavity 
reduces  the  pressure  fluctuations  significantly  inside  the  cavity  for  the  non-reactive  flows.  We  obtained  a  similar 
result  for  the  reactive  flows  with  the  ignition  of  the  fuel  supplied  initially  in  the  cavity. 

The  Injector-cavity  system  was  also  considered  with  the  hydrogen  fuel  of  M  =  1,  constantly  injected  ahead 
of  cavity,  and  we  investigated  how  cavity  plays  a  role  as  a  flame-holder  and  how  the  recessivity  of  cavity  could 
affect  the  flame-holding  efficiency  of  this  system. 

The  results  show  that  the  pressure  fluctuations  in  cavities  with  lower  angle  of  the  aft  are  weaker  than  in 
cavities  with  higher  angles.  It  is  also  shown  that  the  attenuation  of  the  pressure  fluctuations  are  obtained  both 
at  the  center  and  the  middle  of  the  floor  of  the  cavity.  We  also  show  that  the  large  recirculation  2one(s)  formed 
inside  the  cavity  obtained  by  the  lower  order  numerical  scheme  is  induced  not  physically  but  rather  artificially 
due  to  the  heavy  numerical  dissipations.  This  result  shows  that  for  these  sensitive  problems,  high  order  accuracy 
should  be  used  in  order  to  minimize  the  effect  of  the  numerical  dissipation. 

We  also  considered  the  case  of  the  reactive  flows  for  the  90  and  30  degree  aft  walls.  Similar  features  of  the 
pressure  fluctuations  are  shown  as  in  the  non-reactive  flows.  However  the  pressure  fluctuations  are  much  more 
attenuated  for  both  the  90  and  30  degree  walls  than  in  the  non-reactive  cold  flows.  These  results  demonstrate 
that  simulations  of  cold  flows  do  not  necessarily  shed  light  on  the  behavior  of  reactive  flows.  The  shear  layer 
is  becoming  weaker  as  the  degree  of  angle  of  the  aft  wall  is  getting  low  and  the  flow  fields  are  becoming  more 
regularized  for  the  case  of  the  lower  angled  wall. 
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t  =  0.175ms 


t  =  0.275ms 


t  =  0.945ms 


t  =  2.26ms 

Figure  1:  The  water  contour  of  the  reactive  flows:  the  water  density  contours  are  given  in  the  left  figures  for 
the  90°  wall  30°  wall  in  the  right.  From  top  to  bottom  the  instant  times  t  are  0.175ms,  0.275ms,  0.945ms 
and  2.26ms  respectively.  The  maximum  and  minimum  contour  levels  are  0.01  and  0.23  respectively  with  the 
number  of  levels  50. 


Figure  1  shows  the  water  contour  inside  the  cavity  for  the  different  angles  at  different  time.  Here  we  define 
the  region  where  the  flames  are  generated  to  be  same  as  the  region  where  the  water  is  produced.  As  the  Hydrogen 
fuel  is  consumed,  the  water  is  produced  and  starts  to  be  expelled  from  the  cavity  to  the  main  channel.  Figure 
1  shows  that  the  lower  angled  wall  cavity  (30°  in  this  case)  maintains  more  water  than  the  90°  wall  cavity  at 
a  given  time.  The  figure  also  shows  that  the  lower  angled  aft  wall  holds  the  flame  (water  in  this  case)  longer 
than  the  90°  wall  cavity. 

These  results  imply  that  the  flame-holding  efficiency  can  be  increased  by  lowering  the  angle  of  the  aft  wall 
of  the  cavity. 


4  WENO  methods 


These  are  methods  using  either  locally  adaptive,  nonlinear  weighted  stencils  in  performing  high  order  inter¬ 
polations  or  reconstructions,  for  problems  containing  both  shocks  and  complicated  smooth  structures  in  their 
solutions.  Weighted  ENO  schemes  have  been  developed  in  the  last  years  and  applied  successfully  to  simu¬ 
late  shock  turbulence  interactions  to  the  direct  simulation  of  turbulence  ;  to  detonation  shock  dynamics  ,  to 
dynamical  response  of  a  stellar  atmosphere  to  pressure  perturbations  ;  to  overheating  problem  in  multimate¬ 
rial  compressible  flows  ;  to  large-eddy  simulations  ;  to  relativistic  hydrodynamics  equations  ;  to  shock  vortex 
interactions  and  other  gas  dynamics  problems  ,to  incompressible  flow  problems  to  semi-conductor  device  simu¬ 
lation  ;  to  underwater  blast-wave  focusing  ;  to  the  composite  schemes  and  shallow  water  equations  ,  to  real  gas 
computations  ;  to  viscoelasticity  with  fading  memory  ;  to  subpixel  interpolation  in  image  processing  etc. 

In  the  following  we  will  briefly  review  recent  advances  concerning  WENO  (w-eighted  essentially  non-oscillatory) 
methods,  in  [8],  Balsara  and  Shu  have  investigated  very  high  order  finite  difference  WENO  schemes.  The  class 
of  methods  includes  seventh,  ninth,  eleventh  and  thirteenth  order  methods,  stabilized  by  a  suitable  monotonicity 
preserving  bounds.  It  is  verified  that  these  methods  are  indeed  uniformly  accurate  of  the  claimed  high  order 
accuracy,  and  they  give  non-oscillatory  shock  transitions  in  the  presence  of  discontinuities.  These  methods 
have  good  potential  for  situations  where  extremely  high  order  of  accuracy  is  important,  such  as  direct  numer¬ 
ical  simulation  of  turbulence.  Further  investigation  to  analyze  and  improve  these  methods  is  proposed  in  this 
proposal. 

Also  on  WENO  methods,  Hu  and  Shu  [59]  have  constructed  high  order  weighted  essentially  non-oscillatory 
(WENO)  schemes  on  two  dimensional  unstructured  meshes  (triangles)  in  the  finite  volume  formulation.  Third 
and  fourth  order  schemes  are  constructed.  Numerical  examples  are  shown  to  demonstrate  the  accuracies  and 
robustness  of  the  methods  for  shock  calculations.  Further  investigation  to  analyze  and  improve  these  methods 
is  proposed  in  this  proposal. 

In  [78],  Montarnal  and  Shu  have  used  a  recently  developed  energy  relaxation  theory  by  Coquel  and  Perthame 
and  high  order  weighted  essentially  non-oscillatory  (WENO)  schemes  to  simulate  the  Euler  equations  of -real 
gas.  A  relaxation  process  is  performed  for  each  time  step  to  ensure  that  the  original  pressure  law  is  satisfied. 
The  necessary  characteristic  decomposition  for  the  high  order  W’ENO  schemes  is  performed  on  the  characteristic 
fields  based  on  the  simple  7-law,  hence  dramatically  reduces  the  computational  cost.  The  algorithm  only  calls 
for  the  original  pressure  law  once  per  grid  point  per  time  step,  without  the  need  to  compute  its  derivatives  or  any 
Riemann  solvers.  Both  one  and  two  dimensional  numerical  examples  are  shown  to  illustrate  the  effectiveness  of 
this  approach. 

In  [82]  and  [83],  Shu  has  given  two  lectures  notes  to  describe  the  construction,  analysis,  and  application  of 
ENO  (Essentially  Non-Oscillatory)  and  WENO  (Weighted  Essentially  Non-Oscillatory)  schemes  for  hyperbolic 
conservation  laws  and  related  Hamilton-Jacobi  equations. 

In  [97],  Zhou,  Li  and  Shu  first  review  WENO  finite  volume  and  discontinuous  Galerkin  finite  element 
methods,  pointing  out  their  similarities  and  differences  in  algorithm  formulation,  theoretical  properties,  imple¬ 
mentation  issues,  applicability,  and  relative  advantages.  They  then  present  some  quantitative  comparisons  of 
the  third  order  finite  volume  WENO  methods  and  discontinuous  Galerkin  methods  for  a  series  of  test  problems 
to  assess  their  relative  merits  in  accuracy  and  CPU  timing. 

WENO  (Weighted  Essentially  Non-Oscillatory)  finite  difference  methods  are  methods  using  either  locally 
adaptive,  nonlinear  weighted  stencils  in  performing  high  order  interpolations  or  reconstructions,  for  problems 
containing  both  shocks  and  complicated  smooth  structures  in  their  solutions.  The  methods  are  based  on  ENO 
schemes  started  with  the  pioneering  work  of  Harten.  Engquist,  Osher  and  Chakravarthy  [51].  ENO  schemes 
bcised  on  point  values  and  TVD  Runge-Kutta  time  discretizations,  which  are  efficient  for  multi-dimensional 
calculations,  w'ere  designed  in  [84,  85].  Weighted  ENO  schemes  have  been  developed  in  [72]  for  the  third 
order  finite  volume  version  in  ID,  and  in  [64]  for  third  and  fifth  order  finite  difference  version  in  ID  and  multi 
space  dimensions,  with  a  general  framework  for  the  design  of  the  smoothness  indicators  and  nonlinear  weights. 
Later,  second,  third  and  fourth  order  finite  volume  WENO  schemes  for  2D  general  triangulations  have  been 
developed  in  [27]  and  [59].  Very  high  order  finite  difference  WENO  schemes  (for  orders  between  7  and  13) 
have  been  developed  in  [8].  Compact  central  WENO  schemes  have  been  developed  in  [68].  ENO  and  WENO 
have  been  successfully  used  to  simulate  shock  turbulence  interactions  [85,  86,  5];  to  the  direct  simulation  of 
turbulence  [86,  95];  to  detonation  shock  dynamics  [6],  to  dynamical  response  of  a  stellar  atmosphere  to  pressure 
perturbations  [17];  to  overheating  problem  in  multimaterial  compressible  flows  [25];  to  large-eddy  simulations 
[32];  to  relativistic  hydrodynamics  equations  [18]:  to  shock  vortex  interactions  and  other  gas  dynamics  problems 
[23.  48,  49,  61];  to  incompressible  flow  problems  [21,  50,  96]:  to  semi-conductor  device  simulation  [24,  62] 
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;to  Hamilton-Jacobi  equations  [79,  63];  to  magnetohydrodynamics  [65]:  to  wave- vortex  interaction  m  rotating 
shallow  w-ater  [66];  to  the  solution  of  Vlasov-Poisson  equations  [67];  to  underwater  blast- wave  focusing  [69];  to 
the  composite  schemes  and  shallow  water  equations  [70,  71],  to  real  gas  computations  [78];  to  viscoelasticiu 
with  fading  memory  [87];  to  subpixel  interpolation  in  image  processing  [88];  etc.  For  a  short  summary  of  the 
development  and  applications  of  ENO  and  WENO  schemes,  see  [81].  For  a  survey  of  ENO  and  WENO  methods, 
see  the  lecture  notes  of  Shu  [82,  83]. 


Figure  3;  Rayleigh-Taylor,  9th  WENO  Aa:  —  A?/  —  2001  400 ’  soo’  leoo 

As  an  example,  we  can  see  from  Fig.  2  and  Fig.  3  that  a  ninth  order  WENO  scheme  computes  the  Rayleigh- 
Taylor  instability  problem  much  better  than  a  fifth  order  WENO  scheme,  at  about  twice  the  CPU  cost,  on 
the  same  mesh.  This  indicates  the  advantage  of  going  to  high  order  accuracy  for  such  problems.  Such  type  of 
problems  will  be  further  investigated. 
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5  Multi-Domain  Methods 


The  original  motivation  for  the  introduction  of  multi-domain  methods  can  be  found  in  the  restrictions  that  the 
fixed  grids,  required  to  ensure  the  high  spatial  accuracy,  impose.  This  fixed  grid  makes  it  difficult  to  utilize 
adaptivity  and,  for  multi-dimensional  problems,  to  address  problems  in  complex  geometries.  Moreover,  the  use 
of  global  spectral  expansions  makes  it  difficult  to  achieve  a  high  parallel  efficiency  on  contemporary  parallel 
computers. 

Many  of  these  concerns  can  be  overcome  if  one  splits  the  computational  domain  into  a  number  of  geomet¬ 
rically  simple  building  blocks,  e.g.,  squares  and  cubes,  and  then  employs  tensor-product  forms  of  the  simple 
one-dimensional  approximations  as  the  basis  of  an  element  by  element  approximation.  While  this  technique 
opens  up  for  the  use  of  a  highly  non-uniform  resolution  and  the  ability  to  model  problems  in  geometrically 
complex  domains,  it  also  introduces  the  need  to  connect  the  many  local  solutions  in  an  accurate,  stable,  and 
efficient  manner  to  reconstruct  the  global  solution. 

5.1  Patching  Techniques 

The  patching  of  the  local  solutions  in  a  way  consistent  with  the  nature  of  the  hyperbolic  problem  can  be 
performed  in  at  least  two  different  yet  related  ways.  We  shall  refer  to  these  two  different  methods  as  the 
differential  and  the  correctional  method,  respectively. 

To  expose  the  differences  between  the  two  methods,  let  us  consider  the  two  domain  scalar  problem 


du  df(u) 
dt  dx 
dv  df{v) 
dt  ^  dx 


=  0  , 


=  0  , 


X  6  [-1,0] 

X  e  [0, 1]  . 


(5.1) 


To  recover  the  global  solution  U  =  [it,u]  under  the  constraint  that  u{0,t)  =  u(0,t),  the  central  issue  is  how  one 
decides  which  of  the  multiple  solutions  at  i  =  0  takes  preference  and  hence  determines  the  evolution  of  u(0,f) 
and  u(0,  t). 

Provided  that  the  initial  conditions  are  consistent  with  the  continuity  condition  it  will  clearly  remain  contin¬ 
uous  if  we  ensure  that  U((0,  t)  =  U({0,  t).  This  approach,  known  as  the  differential  method,  involves  the  exchange 
of  information  between  the  two  domains  to  ensure  that  the  flux  of  u{0,t)  and  u(0,<)  are  identical  throughout 
the  computation.  There  are,  however,  several  ways  to  do  so. 

In  general  one  introduces  the  flux-derivative 


A  = 


and  requires  that  u  and  v  be  updated  at  x 


du 
0  as 


1/(0, <) 


dv 


vm) 


du 

dt 


x=0 


dv 

dt 


=  -l(A-f-|A|) 


x=0 


du 

dx 


(A-|A|) 


x=0 


dv 

dx 


I=:0 


This  can  be  recognized  as  nothing  else  than  pure  upwinding.  The  extension  to  systems  of  equations  employs 
the  characteristic  form  of  the  system  and  the  multi-dimensional  case  is  treated  by  dimensional  splitting. 

An  alternative  formulation,  based  on  the  weakly  imposed  boundary  conditions  introduced  in  [52,  58,  54], 
takes  the  form 


du 

dt 

dv 

dt 


df{u) 


x  =  0 


-t- 


dt 

df{v) 


dt 


x=0 


r=0 


=  (u(0,  t)  -  i;(0,  t))  , 


which  again  amounts  to  upwinding,  although  on  a  weak  form.  The  advantage  of  this  latter  formulation  is  that 
it  allows  for  establishing  stability  and  it  makes  the  enforcement  of  very  complex  interface  conditions  simple. 
The  extension  to  systems  employs  the  characteristic  variables  and  is  discussed  in  detail  in  [52,  54]  while  the 
multi-dimensional  case  is  treated  in  [58].  Similar  developments  for  methods  employing  multi- variate  polynomials 
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[55,  56]  or  a  purely  modal  basis  defined  on  triangles  and  tetrahedra  has  recently  been  developed,  pawing  the 
way  for  the  formulation  of  stable  spectral  methods  for  the  solution  of  hyperbolic  conservation  laws  using  a  fully 
unstructured  grid. 

Rather  than  correcting  the  local  temporal  derivative  to  ensure  continuity  of  the  flux  across  the  interface 
one  could  choose  to  modify  the  solution  itself.  This  observation  provides  the  basic  foundation  for  correctional 
methods  in  which  both  u  and  v  is  advanced  everywhere  within  the  each  domain,  leading  to  a  multiplicity  of 
solutions  at  x  =  0.  For  the  specific  case  discussed  here,  the  correctional  approach  amounts  to 


u{0,t)  =  v{0,t) 


_  /  u(0, 
1  »^(0, 


t)  if  A  >  0 
t)  if  A  <  0 


which  we  again  recognize  as  upwinding.  The  system  case  is  treated  similarly  by  exploiting  the  characteristic 
variables.  As  for  the  differential  methods,  the  use  of  the  characteristics  implicitly  assumes  a  minimum  degree 
of  smoothness  of  the  solution.  However,  as  no  information  about  the  spatial  derivatives  are  passed  between 
domains,  the  correctional  method  imposes  no  constraints  on  the  smoothness  of  the  grid. 

The  main  appeal  of  the  correctional  method  is  its  simplicity  and  robustness  and  it  has  been  utilized  to 
formulate  very  general  multi-domain  method  for  problems  in  gas-dynamics,  and  in  acoustics  and  elasticity,  and 
electromagnetic. 

Note  that  both  methods  employ  the  local  flux-Jacobian,  A,  which  implicitly  requires  a  certain  amount  of 
smoothness  of  the  solution  at  the  interface.  A  differential  method  overcoming  this  can  be  realized  by  borrowing 
a  few  ideas  from  classical  finite  volume  methods. 

Consider  the  cell  averaged  formulation 


dUj  f{u{xj+i/2))  -  /(u(Xj-i/2)) 

dt  Axj 

where 


_  1 

Axj  =  ii+i/2  -  x,_i/2  ,  =  /  u{s)ds  . 

Here  Xj±i/2  signifies  the  Chebyshev-Gauss-Lobatto  grid  and  xj  refers  to  the  interlaced  Chebyshev-Gauss  grid. 
No  assumptions  are  made  about  the  smoothness  of  the  flux  and  since  each  individual  cell  requires  reconstruction, 
the  patching  of  the  subdomains  is  achieved  by  flux-splitting  techniques  known  from  finite  volume  methods.  This 
approach  was  first  proposed  for  Fourier  methods  and  subsequently  for  the  Chebyshev  approximation  and  has 
the  advantage  of  being  conservative  by  construction.  The  averaging  and  reconstruction  procedure,  which  can  be 
done  in  an  essentially  non-oscillatory  way,  is  essential  for  the  accuracy  and  stability  of  the  scheme  and  several 
alternatives,  exploiting  a  similar  framework,  has  been  also  proposed. 

The  use  of  a  staggered  grid,  collocating  the  solution  u  at  the  Gauss  grid  and  the  fluxes,  f{u),  at  the  Gauss- 
Lobatto  grid,  has  the  additional  advantage  of  allowing  for  the  formulation  of  multi-dimensional  multi-domain 
methods  with  no  grid  points  on  the  vertices  of  the  elements.  This  approach  has  been  developed  for  smooth 
problems  and  eliminates  complications  associated  with  the  treatment  of  vertices  in  multi-domain  methods. 

5.2  Conservation  Properties  of  Multi-Domain  Schemes 

The  important  question  of  the  conservation  properties  of  multi-domain  schemes  is  discussed  in  [14]  in  which  the 
following  polynomial  approximation  to  Eq.(5.1)  is  considered 

^  =  r,Q+{x)  [f+{uN{0,t))  -  /+(u/v(0,t))] 

+T2Q'I'{x)  [/'(u/v(0,«))  -  /■'(■i;;v(0,t))]  -f  SV(un) 

^  =  -rsQjAx)  -  /+(uyv(0,0)] 

-t4Q7;(x)  [/“(uA'(0,t))  - /“(u/v(0,<))]  +  SV{vr^)  , 

where  are  polynomials  that  vanish  at  all  collocation  points  except  x  =  ±1.  Furthermore,  SV{ui\t)  and 

SV{vn)  represent  the  vanishing  viscosity  terms,  or  filtering,  required  to  stabilize  the  nonlinear  problem  as 
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discussed  in  Sec.  2.2.1  and  Sec.  2.2.2,  while  /  =  /+  +  /"  signifies  a  splitting  into  the  upwind  and  downwind 
components  of  the  flux. 

To  establish  conservation  of  the  approximation,  consider  a  test  function  tpix)  and  denote  by  ■!/>/  and  xpji 
its  restriction  to  the  first  and  second  domain  respectively.  We  can  assume  that  tpi  and  iIjji  are  polynomials  of 
order  N  -1  and  that  ipj  vanishes  at  x  =  -1  of  the  first  domain  while  vanishes  x  =  1  of  the  second  domain, 
but  not  at  X  =  0. 

Repeated  integration  by  parts  using  the  fact  that  vanishes  on  the  left  boundary  and  that  SV{un)  vanishes 
at  the  boundaries  of  each  domain  yields 


iPj{x)SV{un)  dx  = 


which  tends  to  zero  with  increasing  N.  A  similar  result  can  be  obtained  for  the  second  domain. 
Consider  now 


To  recover  that  and  vn  are  weak  solutions  to  Eq.(5.1),  i.e.,  the  above  integral  vanishes,  we  must  require 

Ti  +  Ts  =  1  ,  T2  +  T4  =  1  .  (fi-2) 


However,  for  linear  stability  one  can  show  that 

1  1 

ri  >  2  ,  r2  <  -  , 

1  1 

T3  <  2  -  ^4  >  2  , 

are  necessary  and  sufficient  to  guarantee  stability. 

This  leaves  us  with  a  set  of  conditions  under  which  to  design  stable  and  conservative  scheme.  In  particular, 
if  we  choose  to  do  pure  upwinding  at  the  interfaces  by  specifying 


n  =  r4  =  1  ,  To  =  T3  =  0  , 

we  essentially  recover  the  discontinuous  Galerkin  method. 

1 

n  =T2  =  T3  =  Ti  =  -  , 

which  yields  a  marginally  stable  and  conservative  scheme  on  the  form 


^  +  =  -^Qji{x)[f{vN{0,t))  -  f{UN{0,t))]  +  SV{vn)  , 

i.e.,  the  interface  boundary  conditions  are  imposed  on  the  fluxes  /  rather  than  on  the  split  fluxes  and  /~. 

5.3  Computational  Efficiency 

An  interesting  question  pertaining  to  the  use  of  multi-domain  methods  is  how  one  decides  how  many  elements 
and  what  resolution  to  use  within  each  element.  In  pragmatic  terms,  what  we  are  interested  in  is  to  identify 
the  optimal  combination  of  the  order  of  the  polynomial,  N,  and  the  number  of  elements,  K,  needed  to  solve  a 
particular  problem  to  within  a  maximum  error  using  minimum  computational  resources.  On  one  hand,  it  is  the 
high  order  of  the  interpolation  polynomial  that  yields  the  very  accurate  approximation.  On  the  other  hand,  the 
computation  of  derivatives  generally  scales  as  0{N^)  while  the  total  work  scales  only  linearly  with  the  number 
of  elements.  To  develop  guidelines  for  choosing  the  optimal  N  and  K,  consider  a  one-dimensional  wave  problem 
with  a  smooth  solution.  Assume  that  the  approximation  error,  E{N,K),  scales  as 
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where  k  is  the  maximum  wavenumber  in  the  solution,  i.e.,  it  is  proportional  to  the  inverse  of  the  minimum  spatial 
length  scale.  We  shall  require  the  maximum  absolute  error,  E,  is  bounded  as  E  <  exp(— 7).  and  estimate  the 
computational  work  as 

W{N,  K)  =  c,  KN"^  +  C2KN  , 

where  C\  and  C2  are  problem  specific  constants.  Minimizing  the  work  subject  to  the  error  constraint  yields  the 
optimal  values 

7rk 

A  opt  =  7  >  ATopt  =  -rj —  exp 

^''opt 

One  observes  that  high  accuracy,  i.e.,  7  large,  should  be  achieved  by  using  a  large  number  of  modes,  N,  and  not, 
as  one  could  expect,  by  employing  many  subdomains  each  with  a  low  number  of  modes.  For  very  smooth  and 
regular  functions,  where  k  is  small,  or  if  only  moderate  accuracy  is  required  a  multi-domain  approach  may  not  be 
the  optimal  method  of  choice.  On  the  other  hand,  if  the  function  exhibits  strongly  localized  phenomena,  i.e.,  k  is 
large,  one  should  introduce  several  domains  to  minimize  the  computational  burden.  While  these  arguments  are 
loose,  they  indicate  that  an  optimal  choice  of  N  and  K  for  most  problems  seems  to  be  a  few  larger  subdomains, 
each  with  a  reasonable  number  of  modes  to  maintain  an  acceptable  spatial  accuracy. 

These  results  have  been  confirmed  in  computational  experiments  in  [52,  58,  57]  indicating  that  =  8  —  16 
is  reasonable  for  two-dimensional  problems  and  =  4  —  8  is  reasonable  for  three-dimensional  problems.  If  this 
results  in  insufficient  resolution  one  should  generally  increase  the  number  of  domains  rather  than  the  resolution. 
Similar  conclusions  have  been  reached  for  the  analysis  of  spectral  multi-domain  methods  in  a  parallel  setting 
[26]. 
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6  Richtmyer  Meskhov  Instability  (RMI) 

Inertial  Confinement  Fusion  Program  (ICF)  uses  high  energy-  pulse  sources  such  as  X-rays  and  lasers  to  illuminate 
the  target  sphere  in  order  to  achieve  auto  fusion  ignition  and  such  is  the  plan  for  the  future  on  the  National 
Ignition  Facility  (NIF)  at  Lawrence  Livermore  National  Laboratory  (LLNL).  It  is  crucial  to  achieve  uniform 
compression  of  the  target  sphere  in  order  to  obtain  maximum  efficiency  as  the  impulsively  accelerated  non- 
uniform  sphere  surface  by  a  shock  wave  causes  a  non-uniform  pressure  profile  over  the  sphere.  One  can  observe 
a  similar  phenomenon  in  the  mixing  of  fuel  with  oxidants  in  a  SCRAMJET  engine,  the  fuel  jets  are  under 
impulsive  acceleration  by  shockwave  .  The  fuel  mixing  efficiency  is  enhanced  by  the  stretching  of  the  fuel- 
oxidant  interface  and  the  breakup  of  fuels  into  finer  droplets.  The  source  of  these  two  phenomena  are  known 
as  the  Richtmyer-Meshkov  Instability  (RMI)  and  they  are  related  to  the  well  known  fluid  instability  studied 
theoretically  by  Richtmyer  and  experimental!}^  by  Meshkov. 

In  a  word,  the  RMI  results  from  an  impulsively  accelerated  interface  of  materials  with  different  densities 
under  perturbation.  This  form  of  instability  is  different  from  the  closely  related  fluid  instability  known  as  the 
Rayleigh-Taylor  instability  in  which  the  material  interface  is  under  a  constant  acceleration  force  such  as  gravity. 
The  vorticity  generated  by  the  cross  product  of  the  pressure  gradient  and  the  density  gradient  deforms  and 
amplifies  the  interface  perturbation  and  grows  in  time.  The  penetration  of  the  heavier  fluid  into  the  lighter 
fluid  forms  spikes  and  bubbles  and  vice  versa.  A  turbulence  mixing  zone  of  the  two  fluids  is  created  and  grows 
with  time. 

In  order  to  capture  the  shock-interface  interaction  and  the  fine  scale  structures  within  the  turbulence  mixing 
zone,  high  order  methods  are  highly  desirable.  From  the  point  of  view  of  the  numerical  calculation,  we  can  break 
the  RMI  problem  into  two  parts.  First,  we  have  the  issue  of  reliably  calculating  the  motion  of  a  possibly  very 
strong  shock  wave,  and  second  we  have  the  issue  of  reliably  calculating  the  mix  that  ensues  after  this  shock  wave 
accelerates  the  interface.  It  is  in  this  second  area  of  calculating  the  ensuing  mix  where  high  order  numerical 
schemes  offer  unparalleled  efficiency.  This  efficiency  comes  from  the  very  fundamental  fact  that  the  truncation 
error  in  the  differentiation  operators  can  be  made  much  smaller  by  increasing  the  order  of  the  scheme  than  by 
increasing  the  number  of  grid  points.  One  of  the  consequences  of  this  ability  to  rapidly  reduce  the  truncation 
error  is  that  long-time  integration  is  computationally  much  less  expensive  with  high  order  schemes,  than  with 
low  order  schemes.  Further,  when  the  flow  contains  small  scale  features  such  as  in  turbulent  flows,  high  order 
numerical  operators  are  far  more  efficient  than  low  order  operators,  see  [98],  and  adapting  a  low  order  scheme 
can  not  overcome  the  efficiency  of  a  high  order  scheme.  Among  the  high  orders  schemes  considered,  Spectral 
methods  (Spectral)  and  Weighted  Essentially  Non-Oscillatory  finite  difference  schemes  (WENO)  are  considered 
in  this  study.  High  Order  compact  schemes  are  another  candidate  but  was  not  considered  in  this  study.  High 
order,  in  the  sense,  means  an  order  of  accuracy  higher  than  two.  The  numerical  schemes  developed  here  will 
served  as  a  basis  for  the  validation  of  the  results  obtained  from  both  experiments  and  other  numerical  methods. 

The  governing  equations  are  the  two/three-dimensional  Euler  equations  in  Cartesian  coordinates  describing 
the  conservation  of  mass,  momentum  and  energy.  The  physical  domain  is  a  rectangular  domain  with  0  <  x  <  Lx 
and  0  <  y  <  X,  where  Lx  is  the  user  specified  domain  length  in  x  and  A  is  the  wave  length  of  a  single  mode 
perturbation  along  the  interface  separating  two  different  gases  in  the  y  direction.  In  this  study,  the  gases  are 
Xenon  (Xe)  and  Argon  (Ar).  The  specific  heat  ratio  7  is  assumed  to  be  the  same  for  both  gases,  7  =  §• 

The  initial  conditions  (see  figure  4)  is  an  incident  shock  of  Mach  number  M  =  4.46  located  at  Xs  =  0.05cm 
that  travels  downstream  toward  the  interface.  Given  the  shock  Mach  number  AI,  temperature  T  in  the  pre-shock 
region,  and  the  density  of  the  Xenon  gas  pxe-  the  initial  condition  of  the  flow  satisfies  the  Hugoniot-Rankine 
Condition  for  normal  steady  shock,  i.e., 

T2  =T  ,P2=  pXe  7  P2  =  RP2T2  ,  C2  =  }/ 'JP2I P2  7  U2  —  Af  (72  ,  V2  =  0 
Pi  =  P2— -  7P1  =  P2  2  =  U2P2/P1  =  0 

The  subscripts  1  and  2  denote  the  pre-shock  and  post-shock  condition,  respectively.  To  specify  condition  for 
the  moving  shock,  the  shock  speed  s  =  MC2  is  subtracted  from  the  pre-  and  post-  shock  velocity  Ui  and  U2, 
respectively.  Using  the  cgm  units,  the  constant  R  =  Ro/AIxe,  where  Pq  =  8.31441  x  lO’^  ^  is  the  universal 
gas  constant,  Mxe  —  131.29  and  M^r  =  39.948  are  the  molecular  weight  of  Xenon  and  Argon, 

respectively. 

In  the  pre-shock  region,  the  temperature  T  =  296.0  K  and  the  density  of  Argon  and  Xenon  gases  are 
pAr  =  0.89  X  10“^  and  pxe  =  2.9  x  10~®  respectively.  The  pressure  is  assumed  to  be  half  of  the  normal 
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atmospheric  pressure. 

Once  the  pre-  and  post-  shock  states  of  Xenon  gas  are  specified,  the  region  of  Argon  gas  will  be  superimposed 
onto  the  pre-shock  region  replacing  the  Xenon  gas.  The  diffused  interface  between  the  two  gases  is  further 
perturbed  to  form  a  sinusoidal  wave  with  some  finite  thickness  with  amplitude  a  =  1.0  cm  and  wave  length 
A  =  3.6  cm. 


Figure  4:  Initial  Condition  of  the  Richtmyer-Meshkov  instability  problem. 


Solution  Procedure  : 

1.  Periodical  Domain  is  specified  in  the  y  direction  and  symmetry  property  of  the  problem  is  exploited  to 
reduce  the  amount  of  computational  operations  by  half. 

2.  A  combined  Chebyshev  and  Fourier  collocation  methods  is  used  to  discretize  the  Euler  equations  in  space 
yielding  a  Spectral  scheme.  A  WENO-LF-5  fifth  order  finite  difference  scheme  is  also  employed  for  the 
solution  of  the  problem. 

3.  The  high-performance  parallel  library  PseudoPack,  written  by  Bruno  Costa  and  Wai  Sun  Don,  is  utilized 
for  both  Spectral  scheme. 

4.  The  third  order  TVD  Runge  Kutta  method  by  Shu  and  Osher  is  used  to  advance  the  solution,  in  time. 

5.  For  the  Spectral  scheme,  a  lO’th  and  9’th  order  exponential  filter  is  used  for  the  differentiation  and  solution 
smoothing  respectively,  at  each  Runge  Kutta  step  unless  otherwise  specified. 

As  evidenced  from  the  results  of  the  Spectral  and  the  WENO  calculations  shown  below,  the  following  major 
features  of  the  Richtmyer-Meshkov  instability  can  be  observed  (see  figure  5)  at  time  t  =  50  x  10“®s,  namely, 

•  Wave  generated  by  the  shock  refraction  behind  the  gas  interface  in  Box  1. 

•  The  penetration  of  the  heavy  (Xe)  to  light  (Ar)  fluid  causes  the  deformation  of  the  interface  into  a  large 
mushroom  shape  structures  in  Box  2  and  the  opposite  in  Box  5.  They  are  refered  as  Spike  and  Bubble 
respectively,  in  the  literatures.  They  move  in  the  opposite  direction  relative  to  each  other  and  form  a  ever 
larger  turbulence  mixming  zone. 

•  Pressure  wave  along  the  transmitted  shock  in  Box  3. 

•  A  small  jet  and  its  vortical  structure  located  in  Box  4.  The  contact  discontinuity  develops  into  a  more 
complicated  vortical  rollups  in  a  finer  and  long  term  simulation  possibly  caused  by  the  Kelvin-Helmholtz 
instability. 

•  Vortical  rollups  of  the  gaseous  interface  inside  Box  6. 
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Figure  5:  The  numbered  regions  enclose  the  most  prominent  flow  features  of  the  Richmyter-Meshkov  instability 
at  time  t  =  50  x  10“®s. 


The  Mach  number  M,  the  Atwood  number  At  and  the  interface  curvature  play  an  important  role  on  the 
growth  of  perturbed  amplitude  on  the  interface.  In  the  particular  set  of  parameters  studied  here  with  high 
Mach  number  M  =  4.46  and  median  Atwood  number  At  «  0.54,  a  formation  of  triple-shock  configuration  along 
the  interface  indicates  that  shock-interface  interaction  is  in  the  ’’hard”  regime.  A  ’’hard”  regime,  as  quoted 
from  Zaytsev  et  al.  [101]  is  ’’the  propagation  of  secondary  shocks  across  the  flow  that  is  accompanied  by  the 
■formation  of  breaks  and  triple  configurations  on  the  refracted  and  reflected  shocks”.  The  triple-shock  formation 
can  be  observed  easily  in  the  early  time  <«  30  x  10~®s. 


Figure  6;  Density  (Top  Row)  and  V-Velocity  (Bottom  Row)  contour  plot  of  the  Richtmyer-Meshkov  instability 
as  computed  by  the  Spectral  scheme.  Domain  length  in  a:  is  Z-x  =  5  cm.  The  interface  thickness  <5  =  0.2  cm. 
The  final  time  is  t  =  50  x  10“®s.  The  resolution  of  the  Spectral  schemes  are  384x192  (Top  Left),  512x256  (Top 
Right)  and  1024x256  (Bottom  Left). 

The  density  p,  velocity  U  and  V  and  total  energy  E  of  the  solution  of  the  Spectral  scheme  figures  (6)  and 
and  WENO-LF-5  scheme  figures  (7)  and  at  time  t  =  50  x  10“®  s  at  various  resolutions. 

It  can  be  observed  that  the  large  and  median  scale  structures  such  as  transmitted  shock,  shocked-interface 
velocity  and  shock  triple  point  are  basically  in  excellent  agreement  with  each  others.  Some  discrepancies  of 
the  fine  scale  structures  along  the  gaseous  interface,  as  can  be  expected  for  numerical  simulation  of  the  Euler 
equations  which  is  sensitive  to  perturbation  in  nature,  are  observed. 
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Figure  7:  Density  (Top  Row)  and  V-Velocity  (Bottom  Row)  contour  plot  of  the  Richtmyer-Meshkov  instability 
as  computed  by  the  WENO-LF-5scheme.  Domain  length  in  a:  is  =  5  cm.  The  final  time  is  t  =  50  x  10“®s. 
The  resolution  of  the  WENO-LF-5  schemes  are  256x128  (Top  Left),  512x256  (Top  Right)  and  1024x512  (Bottom 
Left). 


The  amplitude  of  the  perturbed  interface  a  is  first  decreased  by  the  compression  of  the  shock  wave.  The 
interface  is  then  accelerated  by  the  shock  and  grows  from  2  cm  to  w  2.6  cm  at  time  143  x  10“®s.  Quantita¬ 
tively,  the  amplitude  growth  rate  a  was  plotted  verus  the  mean  displacement  distance  A'/cia  for  various  initial 
amplutides  Oq  and  Mach  number  M.  It  seems  to  match  the  experimental  data  given  in  figure  1.9  in  the  paper 
by  Zaytsev  et  al.  [101]  for  the  case  of  initial  amplitude  Oq  =  1  cm  and  the  distance  passed  by  shock  after 
interaction  with  the  interface  ~  13  cm.  It  is  consistent  with  the  experimental  observations  for  the  linear 
growth  of  a  within  ’’soft  regular”  regime  and  decrease  growth  of  a  within  the  irregular  regime. 


Figure  8;  Grow'th  rate  of  the  amplitude  (Left)  and  distance  of  interface  traveled  (Right). 
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7  Unstructured  Spectral  and  WENO  Schemes 

While  the  structured  grid  multi-domain  schemes  described  above  are  accurate,  stable  and  robust,  they  suffer 
from  two  main  disadvantages  that  makes  the  use  of  them  difficult  for  certain  classes  of  problems.  It  is  well 
known  that  automated  grid  generation  from  CAD  files  is  very  complex  when  the  geometric  building  blocks  are 
curvilinear  cubes  only.  Moreover,  for  problems  with  strongly  local  behavior,  amenable  to  a  space-time  adaptive 
approach,  a  structured  method  is  complicated  by  the  need  to  propagate  grid  refinements  through  the  grid  or, 
alternatively,  to  use  a  non- conforming  discretizations.  To  directly  address  these  concerns,  we  have  recently 
initiated  a  concentrated  effort  focused  on  the  development  of  methods  for  solving  the  Navier  Stokes  equations 
in  complex  geometries  using  an  unstructured  grid. 

We  built  upon  an  existing  unstructured  grid  high-order  framework,  USEMe  (Unstructured  Spectral  Element 
Method),  in  which  the  computational  problem  is  assumed  described  by  a  collection  of  nonoverlapping  simplices. 
It  is  furthermore  assumed,  as  in  finite  element  methods,  that  the  simplex-based  discretization  is  fully  body- 
conforming  but,  contrary  to  classical  finite  element  methods,  we  allow  the  simplices  to  include  higher  order 
-information  of  the  boundaries  and  interfaces,  i.e.,  the  simplices  are  fully  bodyconforming  with  the  possibility 
of  curved  edges  and  faces.  This  enables  the  use  of  large  higher  order  elements  while  maintaining  a  geometry 
description  to  the  order  of  the  scheme. 

Within  each  of  the  simplices  the  solution  is  assumed  to  be  polynomial  as 

N 

q{x,t)  ^  q^{x,t)  =  'Y^q{xi,i)Li{x)  , 

i=l 

where  Li{x)  represents  the  truly  multidimensional  n’th  order  Lagrange  interpolation  polynomial  based  on  the 
N  =  ^{n  +  l)(n  -f  2)(n  -f  3)  nodal  points,  Xj,  wdthin  the  tetrahedron. 

The  fully  unstructured  grid  approach  with  high-order  representations  of  the  solutions  directly  addresses 
the  need  to  provide  a  scheme  with  high  order  accuracy.  Very  significant  geometric  flexibility  is  furtheripore 
ensured  by  the  possibility  of  utilizing  exiting  geometry  description  and  automated  grid  generation  tools  from 
finite-element  methods. 

It  is  worth  while  emphasizing  that  the  boundary/interface  conditions  are  not  imposed  exactly  but  rather 
weakly  through  the  penalizing  surface  integral  and  that  the  formulation  is  inherently  discontinuous,  enforcing 
the  interface  conditions  weakly  through  the  penalizing  term  and  giving  rise  to  a  highly  parallel  formulation  of 
the  scheme. 

7.1  Weighted  Essentially  Nonoscillatory(WENO)  Methods 

The  other  extreme  is  to  choose  a  local  zeroth  order  basis  in  which  case  qj^  now  takes  the  role  of  the  cell 
average  and  the  high-order  scheme  described  in  the  above  reduces  to  a  classical  unstructured  grid  finite- volume 
scheme.  In  this  scenario  high-order  accuracy  is  achieved  in  the  reconstruction  phase  of  the  fluxes  at  the  element 
interfaces  although  the  particular  choice  of  fluxes  at  the  element  phases  also  severely  impact  the  performance 
of  the  scheme. 

The  central  idea  in  the  weighted  essentially  non-oscillatory(WENO)  methods  is  to  use  locally  adaptive, 
nonlinear  weighted  stencils  to  perform  the  high  order  reconstructions,  borrowing  ideas  from  ENO  schemes. 
Both  ENO  and  WENO  methods  have  proven  themselves  very  successful  in  solving  problems  dominated  by 
strong  shocks  but  with  significant  small  scale  dynamics,  e.g.,  mixing  problems. 
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8  Wellposed  Perfectly  Matched  Layers  for  Advective  Acoustics 


The  ability  to  simulate  accurately  wave  phenomena  is  important  in  several  physical  fields,  e.g.,  electromagnetics, 
ambient  acoustics,  advective  acoustics  associated  with  a  mean  flow,  elasticity  and  seismology. 

Often  the  numerical  simulations  of  such  problems,  due  to  limited  computing  resources,  must  be  confined 
to  truncated  domains  much  smaller  than  the  physical  space  over  which  the  wave  phenomena  take  s  place.  In 
such  cases,  numerical  reflections  of  outgoing  waves  from  the  boundaries  of  the  numerical  domain  can  falsify  the 
computational  results.  This  artifact  limits  the  overall  order  of  accuracy  of  the  algorithm  used  in  the  computation. 
This  is  particularly  troublesome  in  cases  where  higher  order  of  accuracy  is  required  by  mode  resolution,  storage 
availability,  etc. 

To  deal  with  these  type  of  problems  local  non-reflecting  boundary  conditions  were  derived  for  the  wave 
equation  by  Engquist  and  Majda  [22]  and  later  by  Bayliss  and  Turkel  [9].  In  practice,  however,  their  effective¬ 
ness  was  limited.  The  notion  of  perfectly  matched  absorbing  layers  (PML)  was  introduced  in  the  context  of 
computational  electro-magnetics  by  Berenger  [11].  The  idea  behind  PML  is  to  attribute  to  the  layers  “material” 
properties  that  modify  the  original  field  equations  so  that  the  waves  will  decay  in  all  directions  of  propagation 
in  the  layers.  However,  it  has  been  shown  by  Abarbanel  and  Gottlieb  [1]  that  the  splitting  technique  used  by 
Berenger,  results  in  a  set  of  P.D.E’s  which  are  only  weakly  well  posed,  i.e.,  they  may  become  ill  posed  under 
certain  perturbations  -  an  example  of  which  is  provided  in  [1].  We  have  suggested  in  [3]  a  new  approach,  to 
derive  the  PML  equations  from  purely  mathematical  considerations.  This  approach  yields,  in  the  CEM  case,  a 
stable  family  of  PML  equations. 

It  is  straightforward  to  show  that  there  is  a  one  to  one  correspondence  between  the  case  of  two-dimensional 
ambient  acoustics  (no  mean  flow)  and  transverse  electromagnetic  wave  propagation.  Hence  the  PML  procedures 
of  CEM  may  be  applied  to  the  case  of  ambient  acoustics. 

In  [2]  we  considered  the  case  of  advected  acoustics  (non-vanishing  mean  flow)  with  the  PML  layers  being 
normal  and  parallel  to  the  mean  flow,,  we  used  an  extension  of  the  procedure  used  in  [3]  to  construct  the  PML 
equations  for  the  case  of  electromagnetics  by  means  of  a  mathematical  procedure. 

We  considered  the  propagation  of  waves  induced  in  a  uniform  two  dimensional  subsonic  flow,  (uoiO),  of  a 
compressible  fluid,  by  small  perturbations.  This  phenomenon  is  described  by  the  linearized  Euler  equations  for 
the  perturbations  of  the  density,  p' ,  and  velocities,  u'  and  v'  as 


dp'  _  dp'  _  du'  dv'  - 
-^  +  uo~^  +  p,—  +  Po^  =  Q 


(6.1) 


du'  du'  cl 

at  ox  pQ 


dv'  _  dv'  Cq 
at  ox  Pq 


(6.2) 

(6.3) 


where  we  assumed  isentropy  of  the  flow,  i.e  po  =  Po(Po  )•  The  speed  of  sound  cq  is  given  by  Cq  =  dp^/dp^  where 
are  the  unperturbed  pressure  and  density  of  the  flow.  The  dimensional  time  and  distances  are  given  by 
t,x  and  y. 

We  non-dimensionalized  this  set  of  equations  by  using  a  reference  length  X2  =  y2  =  L  (usually  the  relevant 
size  of  the  computational  domain  and  the  wavelength  of  the  illuminating  wave),  and  a  reference  time  tr  =  L/cq. 
Similarly  p^  =  pQ  and  u^  =Vr  cq.  With  M  =  uq/^,  the  resulting  set  of  dimensionless  equation  is: 


Pt  +  Mpx  +Ui  +  Vy  =  0  ,  (6.4) 

Ut  -f  Mwi  -I-  px  =  0  ,  (6.5) 

U; -t- A/Ux +Py  =  0  ,  (6.6) 

where  the  prime  ((•)')  has  been  dropped  from  the  perturbation  quantities.  The  case  of  ambient  acoustics  is 
obtained  by  letting  the  Mach  number  M  — ^  0.  This  case  has  been  discussed  by  Hesthaven  [53],  and  is  known  to 
correspond  exactly  to  the  case  of  2-dimensional  electromagnetic.  Hence,  for  M  =  0,  the  solution  of  any  smooth 
initial  boundary  value  problem  can  be  shown  to  be  a  superposition  of  plane  waves  on  the  form 


(6.7) 


26 


resulting  in  a  dispersion  relation  on  the  form 


a^  +  B-  =  l  .  (6.8) 

When  M  ^  0,  however,  the  resulting  dispersion  relation  is  much  more  complicated,  and  the  analysis  from  the 
case  of  electromagnetics  cannot  easily  be  carried  over  to  the  case  acoustics  waves. 

To  overcome  this  difficulty  equations  (6.4)-(6.6)  were  transformed  to  a  new  set  of  coordinates, 

i  =  x  ,  (6.9) 

Tj  =  y/l  -  M^y  =  72/  ,  (6.10) 

r  =  Mx  +  .  (6.11) 

This  transformation  is  related  to  the  one  utilized  in  [9]  although  with  stretching  applying  in  y  rather  than  in  x 
as  in  [9].  As  we  shall  see  shortly,  this  difference  is  crucial.^ _ 


The  transformed  equations  take  the  form  (with  7  =  \/l  —  M^), 

Vt  +  Mv^  +  ypy  =  0  , 

(6.12) 

M 

(6.13) 

Ut  ^77  —  0  , 

7 

pT  +U^  +  ^Vr,  =  0  . 

(6.14) 

Note  that  the  order  of  the  equations  has  been  reorganized  such  that  that  for  M  =  0,  (7  =  1),  one  recovers 
the  two-dimensional  transverse  electric  set  of  Maxwells  equations  [92]  through  the  simple  transformation,  p 
H,u  Ey,v  -Ex--  This  is  done  purely  for  convenience. 

The  plane  wave  solutions  in  the  stretched  space  (^,  77,  r)  are  of  the  form 

(6.15) 

For  this  ansatz  to  be  a  solution,  A  must  be  the  solvability  eigenvalue  of  (6.12)-(6.14)  after  the  substitution  of 
(6.15).  The  three  distinct  eigenvalues  are  given  as 


Ao  =  1/M 


(6.16) 


Ai  =  Vl  -  =  A  , 

A2  =  -\/l-B2  =  -A  . 


with  the  three  corresponding  eigenvectors 


being 


/  Bj 
A- M 
\  l-MA 


Bl  \ 
-A-M 
1-+-MA  / 


(6.17) 

(6.18) 


(6.19)(a,  6,  c) 


The  Ao-solution  corresponds  to  the  rightward  moving  vorticity  wave  whose  amplitude  tends  to  zero  as  M  — t  0, 
see  eq.  (6.19a).  The  Ai  and  A2  solutions  represent  the  two  counter-propagating  acoustic  waves  moving  to  the 
right  and  left,  respectively,  in  the  (^,72)-plane.  Note  that  because  of  the  specific  transformation  x,y.,t  — >  ^,72,r 
the  eigenvalues  Ai  =  A  =  — A2  satisfy  the  standard  dispersion  relation 


=  1  ,  (6.20) 

analogous  to  (6.8).  Note  also  that  in  the  physical  plane  the  expression  =  1  does  not  constitute  a 

dispersion  relation. 

The  set  of  equations  (6.12)-(6.14)  is  to  be  solved  on  a  finite  computational  domain  contrary  to  the  original 
analytical  problem  which  is  set  on  an  infinite  domain.  We  would  like  to  ensure  that  waves  leaving  the  domain  are 
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not  reflected  back  into  it  as  that  could  otherwise  interact  with  the  solution  and  hence  falsify  it.  The  approach 
taken  here  is  to  surround  the  computational  domain  with  finite  width  strips  which  must  be  defined  such  that 
the  waves  propagate  into  these  absorbing  layers  without  reflection,  and  decay  as  they  continue  their  travels 
inside  these  layers.  Moreover,  we  shall  require  that  these  properties  are  independent  of  the  frequency  as  well  as 
the  angle  of  incidence  of  the  incoming  wave. 

We  suggested  the  following  Layers: 

The  PML  equations  in  the  absorbing  “inflow”  ^-layer  are 


Vr  +  Mv^+  'fPr)  =  -(JxV  +  ‘^<yxQx  +  o'iRx  +  ,  {6-21) 

M 

Ut-  +  Pc  -  — =  -CTxU  ,  (6.22) 

(6.23) 

(6.24) 

(6.25) 

(6.26) 

Note  that  from  a  computational  point  of  view  (6.22)-(6.25)  hardly  add  to  the  amount  of  computing.  The 
quantity  dpi  dp  is  evaluated  in  (6.22)  and  thus  (6.24)-(6.26)  weigh  as  three  additional  o.d.e’s  rather  than  3 
additional  p.d.e’s.  Transforming  the  system  (6.22)-(6.25)  back  to  the  physical  space  (x,y,t)  yields 


Pr  +u^  +  -Vr,  =  -CTxP  , 
7 


-7Pv 


dQx 

dr 


dP.  p 

—  =  V-  OxPx  , 


dR, 

dr 


Qx 


dv  ^^dv  dp  n  ^ 


+  (7^R, 


+  Ma'^Px 


du  du  dp 


dp  ^,dp  du  dv 
dt  dx  dx  dy 


-axP  -  (JxMu  , 


dQx  2 


dPx 

dt 


=  7^(u  -  aPx) 


Let  us  briefly  consider  the  issue  of  wellposedness  of  this  new  set  of  equations.  Clearly,  since  the  equations  for 
Px  and  Rx  are  o.d.e.’s  these  have  no  effect  on  the  issue  of  wellposedness.  The  equation  for  Qx,  however,  may 
affect  the  wellposedness  of  the  original  set  of  equations,  but  we  showed  that  it  did  not. 

The  PML  equations  for  both  the  upper  and  lower  77-layers  of  the  form 


Vr  T  Ajv^  -f-  'yPrj  -  ^OyQy  -|-  OyRy  -f-  'yUyPy  , 


M 

Ur  +  Pc - '"'7 

7 

1 

Pt  +  Me  d - 

7 


=  0  , 
=  0  , 


9Qy 

dr 


=  -  2rQ, 


r^R„ 


d^ 

dr 


{p~aPy)  , 


<^vPy 
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In  the  (a:,  j/,  f)-space  the  system  is 


dr 


=  Qy  ■ 


du  „,du  dp 


9p  ,,dp  du  dv 

—  +  M—  -\ - 1 - 

dt  dx  dx  dy 


0  , 


dQy 

dt 


~  ^  ^  dy  ^y^y 


(6.27) 


dt 


■=  l‘^[p  -  (^yPy] 


Wellposedness  of  this  system  follows  directly  for  the  observation  that  only  spatial  derivatives  of  the  density  is 
introduced  which,  as  we  saw  for  the  system  (6.27)  for  the  ^-layer,  does  not  affect  wellposedness. 

The  development  of  efficient  and  accurate  absorbing  boundary  conditions  for  problems  in  acoustics  and 
beyond  remains  a  very  signihcant  challenge.  What  we  have  presented,  however,  provided  a  mathematical 
framework  in  which  such  development  may  be  successful.  Indeed,  the  development  of  a  PML  for  the  three- 
dimensional  equations  of  acoustics  is  straightforward  provided  only  that  the  mean  flow  can  be  considered 
spatially  constant. 

Of  equal  importance  is  the  development  of  PML  methods  for  problems  involving  smoothly  varying  mean 
flows,  as  in  boundary  layers  and  jets.  While  the  mathematical  tools  developed  so  far  certainly  are  applicable 
for  sufficiently  smooth  variations,  new  developments  are  most  likely  needed  to  address  the  general  variable 
coefficient  problem. 
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